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Q I Abstract. In this study we have computed the pair correlation functions in 

^_J ' the two-dimensional Hubbard model using a quantum Monte Carlo method. We 

employ a new diagonalization algorithm in quantum Monte Carlo method which is 
free from the negative sign problem. We show that the d-wave pairing correlation 
function is indeed enhanced slightly for the positive on-site Coulomb interaction 



C^ 



^ ■ 

^ , U when doping away from the half-filling. When the system size becomes large, 

the pair correlation function Pd is increased for U > compared to the non- 
interacting case, while P^j is suppressed for U > when the system size is small. 
The enhancement ratio Pd[U]/Pd[U = 0] will give a criterion on the existence 
of superconductivity. The ratio Pd[U]/Pd[U = 0] increases almost linearly oc L 
as the system size L X L is increased. This increase is a good indication of 

•^J . an existence of superconducting phase in the two-dimensional Hubbard model. 

C ' There is, however, no enhancement of pair correlation functions at half-filling, 

O ' which indicates the absence of superconductivity without hole doping. 
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1. Introduction 

Strongly correlated electron systems have been studied intensively in relation to high- 
temperature superconductivity (SC). High-temperature superconductors [J El |3l S] 
are known as a typical correlated electron system. Recently, the mechanism of 
superconductivity in high-temperature superconductors has been extensively studied 
using various two-dimensional (2D) models of electronic interactions. Among them 
the 2D Hubbard model[5] is the simplest and most fundamental model. This model 
has been studied intensively using numerical tools, such as the quantum Monte Carlo 
(QMC) method [Sl[71ll[Il[TnilIIllIl[Il[Il[Il[ISl[ni[IHl[ini[2ni[2I], and the 
variational Monte Carlo (VMC) method[ll [SaiMlllillinillZllSiEllSnilSIllSlISS]- 

The Quantum Monte Carlo method is a numerical method employed to simulate 
the behavior of correlated electron systems. It is well known, however, that there 
are significant issues associated with the application of the QMC method. The 
most important one is that the standard Metropolis (or heat bath) algorithm is 
associated with the negative sign problem. In past studies workers have investigated 
the possibility of eliminating the negative sign problemflHl [ITl [T9[|21| . 

In this paper we adopt an optimization scheme which is based on diagonalization 
Quantum Monte Carlo (QMD) methodlH] (a bosonic version was developed in 
Ref . |34) ) ■ as well as the Metropolis Quantum Monte Carlo method (called the 
Metropolis QMC in this paper). In general, and as in this study, the ground-state 
wave function is defined as 

^ - e-^"%, (1) 

where H is the Hamiltonian and tpQ is the initial one-particle state such as the Fermi 
sea. In the QMD method this wave function is written as a linear combination of 
the basis states, generated using the auxiliary field method based on the Hubbard- 
Stratonovich transformation; that is 



^ = E' 



(2) 



where </>„ are basis functions. In this work we have assumed a subspace with 
^states basis wave functions. From the variational principle, the coefficients {c„i} 
are determined from the diagonalization of the Hamiltonian, to obtain the lowest 
energy state in the selected subspace {(f>m}- Once the Cm coefficients are determined, 
the ground-state energy and other quantities are calculated using this wave function. 
If the expectation values are not highly sensitive to the number of basis states, we 
can obtain the correct expectation values using an extrapolation in terms of the basis 
states in the limit N states — ^ oo. 

Whether the 2D Hubbard model can account for high-temperature superconduc- 
tivity is an important question in the study of high-temperature superconductors. 
In correlated electron systems, there is an interesting phenomenological correlation 
between the maximum T^ and the transfer integral t: 

ksTc ~ O.U/{m*/m). (3) 

m*/m indicates the mass enhancement factor and t^ff = t/{m*/m) is the effective 
transfer integral. By adopting t ^ 0.5eV[35] and m* /m ~ 5, this formula applies to 
high- Tc cuprates with Tc ~ lOOK. As the electron becomes heavier, Tc is lowered (in 
accordance with the lowering of Tc in the underdoped region). We can choose t ^ O.leV 
and m * /m ~ 2 for iron pnictides to give T^ ~ 50K. This formula strongly suggests 
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that high-temperature superconductivity originates from the electron correlation, not 
from the electron-phonon interaction. 

Most of QMC method results do not support superconductivity, although 
the results of VMC method with the Gutzwiller ansatz indicates the stable d- 
wave pairing state for large U. The computations of the pair-field susceptibility 
suggest the existence of the Kosterlitz-Thouless transition in the 2D Hubbard model 
indicating superconducting transition in real 3D systems [36l |37]. The perturbative 
and Random phase approximation (RPA) calculations also support superconductivity 
with anisotropic pairing symmetry [5S1 [321 SS SIl EH- In contrast, the pair 
correlation functions obtained by a QMC method [18] are extremely suppressed for the 
intermediate values of U. This result suggests that superconductivity is impossible in 
the 2D Hubbard model. The objective of this paper is to compute pair correlation 
functions and clarify this discrepancy using a new QMC method with employing the 
diagonalization scheme pTj. We show that the pair correlation function is indeed 
enhanced at doping. 

2. Model and the Wave function 

2.1. Hamiltonian 

The Hamiltonian is the Hubbard model containing on-site Coulomb repulsion and is 
written as 

H= ~Y1 *'J (^t^J- + ^■^•) + ^Y1 "J't"ji, (4) 

where c^^ (cja) is the creation (annihilation) operator of an electron with spin a at 

the j-th site and Uja- = clj^Cja. tij is the transfer energy between the sites i and j. 
tij — t for the nearest-neighbor bonds and t^ = ~t' for the next nearest-neighbor 
bonds. For all other cases tij =0. U is the on-site Coulomb energy. The number of 
sites is N and the linear dimension of the system is denoted as L, i.e. N = Li^ . The 
energy unit is given by t and the number of electrons is denoted as N^ . 

2.2. Quantum Monte Carlo method - Metropolis algorithm 

In a Quantum Monte Carlo simulation, the ground state wave function is 

V- = e-^^Vo, (5) 

where "00 is the initial one-particle state represented by a Slater determinant. For 
large r, e""^^ will project out the ground state from tpQ. We write the Hamiltonian as 
H = K -\-V where K and V are the kinetic and interaction terms of the Hamiltonian 
in Eq.Q, respectively. The wave function in Eg.© is written as 

for r = At • m. Using the Hubbard-Stratonovich transformational 143] . we have 

exp(-ATt/ni^ni|.) = - ^ exp(2aSj(ni^-ni^)--t/Ar(nj^+ni^)), (7) 

Si=±l 



(8) 
V„{{s,}) ^ 2aaJ2s^n^a - ::UATJ2n^a, (10) 
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for (tanha)^ — tanh(ATL//4) or cosh(2a) = e'^'^^^^. The wave function is expressed 
as a summation of the one-particle Slater determinants over all the configurations of 
the auxiliary fields Sj — ±1. The exponential operator is expressed as[i5] 

(^-A.K^-A.v)™ ^ _i_ j2 n5™(..M) 

{sum '^ 

xB^_,{s,{m-l))---B1{s,{l)), 
where we have defined 
for 

I I 

Ka^-Y,t,j{clcj^+h.C.). (11) 

ij 

The ground-state wave function is 

■4^ ^^Cn(f>n, (12) 

n 

where 4>n is a Slater determinant corresponding to a configuration {si{£)} {i = 
1, • • • , A'^; ^ = 1, • • • , m) of the auxiliary fields: 

0n=n^™(s.("^))-- -51^(5. (1))^0 

(7 

^ €^l (13) 

The coefficients c„ are constant real numbers: Ci — C2 — ■ ■ ■■ The initial state V'o is 8- 
one-particle state. The matrix of Va{{si}) is a diagonal matrix given as 

y<T{{si}) = diag(2aCTSi - [/At/2, • • • , 2acrsjv - C/At/2). (14) 

The matrix elements of K^ are 

{Kcr)ij = ~ t i,j are nearest neighbors 

= otherwise. (15) 

(f)'^ is an A^ X N^r matrix given by the product of the matrices e"'^'^^", e^" and ^q . 
The inner product is thereby calculated as a determinant [17]. 

(<^^C) = det(^,^t^,-j_ (16) 

The expectation value of the quantity Q is evaluated as 

Pgn = det(0£ 05^)det(0^'^(/)^°') can be regarded as the weighting factor to obtain 
the Monte Carlo samples. Since this quantity is not necessarily positive definite, the 
weighting factor should be |Pf„|; the resulting relationship is, 

{Qa)= 5]^^"(Q^)^«/^Pm 
in in 

= ^ \Pin\sign{Pen){Qa)in/ y^ \Pin\sign{Pen) 

in in 

(18) 
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where sign{a) = a/\a\ and 

This relation can be evaluated using a Monte Carlo procedure if an appropriate 
algorithm, such as the Metropolis or heat bath method, is employed|43j. The 
summation can be evaluated using appropriately defined Monte Carlo samples, 

(0,) = ^S'-'^''"'^"->'Qf- . (20) 

where Umc is the number of samples. The sign problem is an issue if the summation 
of sign(Pen) vanishes within statistical errors. In this case it is indeed impossible to 
obtain definite expectation values. 

2.3. Quantum Monte Carlo method - Diagonalization algorithm 

Quantum Monte Carlo diagonalization (QMD) is a method for the evaluation of (Qa) 
without the negative .sign problem. The configuration space of the probability [|Pmn|| 
in Eq. (120)) is generally very strongly peaked. The sign problem lies in the distribution 
of Pmn in the configuration space. It is important to note that the distribution of 
the basis functions 4>rn ("^ = 1,2,---) is uniform since Cm are constant numbers: 
ci = C2 = ■ ■ ■■ In the subspace {(j^m}, selected from all configurations of auxiliary 
fields, the right-hand side of Eq. (jl7[) can be determined. However, the large number 
of basis states required to obtain accurate expectation values is beyond the current 
storage capacity of computers. Thus we use the variational principle to obtain the 
expectation values. 

From the variational principle, 

/Q\ _ Z-^mn CmCn{(PmQ(Pn) .^IN 

where c™ {m = 1, 2, • • •) are variational parameters. In order to minimize the energy 
the equation dEjdcn = (n := 1, 2, • • •) is solved for, 

i{4>nH(j)yn) - Ey^Cyyi{(f)n4>ra) = 0. (23) 



II c„ 



If we set 

i?m„ = {^rnH(f>n), (24) 

A-mn = {4>m(t>n), (25) 

the eigen equation is 

Hu = EAu, (26) 

for u = (ci, C2, • • •)*. Since (j)„i {m = 1,2,---) are not necessarily orthogonal, A is 
not a diagonal matrix. We diagonalize the Hamiltonian A~^H, and then calculate 
the expectation values of correlation functions with the ground state eigenvector; in 
general A~^H is not a symmetric matrix. 
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Figure 1. Pair correlation function Dyy{£) and Dyx{(-) for 4 X 3, (7 = 4 and 
Ne = 10 obtained by tlie diagonalization quantum Monte Carlo method (a) and 
the Metropolis quantum Monte Carlo method (b). The square are the exact 
results obtained by the exact diagonalization method. In (a) the data fit using a 
straight line using the least-square method as the variance is reduced. We started 
with Nstates = 100 (first solid circles) and then increase up to 2000. 



In order to optimize the wave function we must increase the number of basis states 
{0m}- This can be simply accomphshed through random sampling. For systems of 
small sizes and small U , we can evaluate the expectation values from an extrapolation 
of the basis of randomly generated states. The number of basis states is about 2000 
when the system size is small. For systems 8x8 and 10 x 10, the number of states in 
increased up to about 10000. 

In Quantum Monte Carlo simulations an extrapolation is performed to obtain 
the expectation values for the ground-state wave function. The variance method 
has been proposed in variational and Quantum Monte Carlo simulations, where the 
extrapolation is performed as a function of the variance. An advantage of the variance 
method lies is that linearity is expected in some cases J441[T9]: 

(Q) - Qexact fX V, (27) 

where v denotes the variance defined as 

„. {{H'{H)?) 



and Qexact is the expected exact value of the quantity Q. 

3. Pair correlation functions 

In this section, we present the results obtained by the QMC and QMD methods. 



(28) 



3.1. Comparison of two methods 

The pair correlation function Dap is defined by 



(29) 
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Figure 2. Pair correlation function Dyy{£) as a function of the energy variance 
V in (a) and l/m in (b) for 30 X 2, U = 4 and Ne = 48. We used (a) the 
diagonalization quantum Monte Carlo method and (b) the Metropolis quantum 
Monte Carlo method. We set the open boundary condition. Prom the top, 
e = (1,0), (2,0), (5,0), (4,0), (3,0) and (6,0). 



where Aa{i), a = x,y, denote the annihilation operators of the singlet electron pairs 
for the nearest-neighbor sites: 

Aa(i) = Ci^Ci+at - CjtQ+aJ.- (30) 

Here a is a unit vector in the a{— x, j/)-direction. We consider the correlation function 
of d-wave pairing: 

Prf(^) = (Arf(z + ^)tArf(i)), (31) 

where 

Arf(i) = A,(i) + A_,(i) - l\y{i) - A-a(i). (32) 

i and i + (, denote sites on the lattice. 

We show how the pair correlation function is evaluated in quantum Monte Carlo 
methods. We show the pair correlation functions Dyy and Dyx on the lattice 4x3 
in Fig.l. The boundary condition is open in the 4-site direction and is periodic in 
the other direction. An extrapolation is performed as a function of l/m in the QMC 
method with Metropolis algorithm and as a function of the energy variance v in the 
QMD method with diagonalization. We keep Ar a small constant ~ 0.02 ~ 0.05 and 
and increase r = Ar • m, where m is the division number m of the wave function 
•ip in eq.(5). In the Metropolis QMC method, we calculated averages over 5 x 10^ 
Monte Carlo steps. The exact values were obtained by using the exact diagonalization 
method. Two methods give consistent results as shown in figures. All the Dyy{i) and 
Dyx{t) are suppressed on 4 x 3 as C/ is increased. In general, the pair correlation 
functions are suppressed in small systems. 

In Fig. 2, we show the inter-chain pair correlation function Dyy{£) as a function 
of l/m (b) and the energy variance (a) for the ladder model 30 x 2. We use the 
open boundary condition. The boundary condition is not important for our purpose 
to check the consistency between QMC and QMD mthods. The number of electrons 
is Ne = 48, and the strength of the Coulomb interaction is C/ = 4. Aj,(i) indicates 
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Figure 3. Pair correlation function P^ as a function of the energy variance v 
on 8 X 8 lattice. U = 3, t' = —0.2 and the electron number is A^e = 54. We have 
shown Pd{t-) = {A(;(i + iy A{i)) for £ = {m,n) — i and i = (1, 1), where {m,n) 
are shown in the figure. 
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Figure 4. Pair correlation function Pj as a function of the distance R = \£\ on 
8x8 lattice for (a) the half-filled case Ne = 64 and (b) Ne = 54. We set t' = 0.0 
and U = 0, 3 and 4 for (a) and t' = -0.2 and U = 0, A and 6 for (b). To lift the 
degeneracy of electron configurations at the Fermi energy in the half-filled case, 
we included a small staggered magnetization ~ 10~* in the initial wave function 
^0- 



the electron pair along the rung, and Dyy {£) is the expectation value of the parallel 
movement of the pair along the ladder. The results obtained by two methods are in 
good agreement except £ = (1,0) (nearest-neighbor correlation). 
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Figure 5. Pair correlation function Pj as a function of f/ on 8 X 8 lattice, 
t' = -0.2 for Ne = 54 (diamonds), and t' = for A^e = 50 (squares) and Ne = 64 
(circles). We have shown P^ (^) = (A^(i+£)t A(j)) for £ = (m,,n)-i and i = (1,1), 
where (m, n) are shown in the figure. 
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Figure 6. Pair correlation function P^ as a function of the distance R = \£\ 
on 10 X 10 lattice for A^e = 82 and t' = -0.2. The strength of the Coulomb 
interaction is C/ = 0, 3 and 5. 
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Figure 7. Enhancement ratio of pair correlation function Pd\u/Pd\u=0 as a 
function of tlie linear system size L for f/ = 4 and (7 = 2. The electron density 
Uf, is about 0.8: rie ~ 0.8 for squares. The data for U = 4 and n^ ~ 0.18 are also 
shown by circles. 
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Figure 8. Enhancement ratio of pair correlation function Pd\u /Pd\u=0 as a 
function of the electron density ng. We adopt t' = —0.2 and (7 = 4. For the 
half-filled case, the diamonds show that for t' = on 8 X 8 lattice (solid diamond) 
and 6x6 lattice (open diamond). The square is for t' = —0.2 on 8 X 8 and 10 X 10 
where there is no enhancement. 
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3.2. Pair correlation in 2D Hubbard model 

We present the results for pair correlation in the two-dimensional Hubbard model. In 
this section we show the results using the diagonalization QMC method because the 
Metropolis QMC method has a negative sign problem. We first examine the 8x8 
lattice. The Pd was estimated by an extrapolation as a function of the variance v, as 
shown in Fig. 3, where the computations were carried out on 8 x 8 lattice with [/ = 3, 
t' = —0.2 and Ne = 54. The extrapolation was successfully performed for 8x8. 

We consider the half-filled case with t' — 0; in this case the antiferromagnetic 
correlation is dominant over the superconductive pairing correlation and thus the 
pairing correlation function is suppressed as the Coulomb repulsion U is increased. 
The Fig. 4(a) exhibits the d-wave pairing correlation function P^ on 8 x 8 lattice as a 
function of the distance. The Pd is suppressed due to the on-site Coulomb interaction, 
as expected. Its reduction is, however, not so considerably large compared to previous 
QMC studies [TSj where the pairing correlation is almost annihilated for U — 4. We 
then turn to the case of less than half-filling. We show the results on 8 x 8 with electron 
number N^ = 54. We show Pd as a function of the distance in Fig. 4(b) {Ne = 54). In 
the scale of this figure, Pd ior U > is almost the same as that of the non-interacting 
case, and is enhanced slightly for large U. Our results indicate that the pairing 
correlation is not suppressed and is indeed enhanced by the Coulomb interaction U, 
and its enhancement is very small. The Fig. 5 represents Pd as a function of U for 
iVe = 54, 50 and 64. We set t' = for N^ = 50 and t' = -0.2 for N^ = 54 so that we 
have the closed shell structure in the initial function. In the system of this size, the 
effect of the inclusion oi t' ^ is small. The Fig. 6 shows Pd on 10 x 10 lattice. This 
also indicates that the pairing correlation function is enhanced for U > 0. There is a 
tendency that Pd is easily suppressed as the system size becomes small. We estimated 
the enhancement ratio compared to the non-interacting case Pdii)\u / Pdi^)\u=o at 
\£\ r^ L/2 for He ^ 0.8 as shown in Fig. 7. This ratio increases as the system size is 
increased. To compute the enhancement, we picked the sites, for example on 8 x 8 
lattice, £ = (3,2), (4,0), (4,1), (3,3), (4,2), (4,3), (5,0), (5,1) with |£| ~ 4 - 5 and 
evaluate the mean value. In our computations, the ratio increases almost linearly 
indicating a possibility of superconductivity. This indicates Pd{£) ^ LPd{i) ^ £Pd{(-) 
for l^ L. Because Pd{t)\u=o ~ l/|^|^ we obtain Pd{(.) ~ £Pd(€) ~ l/|fp for \i\ - L. 
This indicatesthat the exponent of the power law is 2. When t/ = 2, the enhancement 
is small and is almost independent of L. In the low density case, the enhancement 
is also suppressed being equal to 1. In Fig. 8, the enhancement ratio is shown as a 
function of the electron density n^ for U — A. A dome structure emerges even in small 
systems. The square in Fig. 8 indicates the result for the half-filled case with t' = —0.2 
on 8 X 8 lattice. This is the open shell case and causes a difficulty in computations as 
a result of the degeneracy due to partially occupied electrons. The inclusion of t' < 
enhances Pd compared to the case with i' = on 8 x 8 lattice. Pd is, however, not 
enhanced over the non-interacting case at half-filling. This also holds for 10 x 10 lattice 
where the enhancement ratio ^^ 1. This indicates the absence of superconductivity at 
half-filling. 

4. Summary 

The quest for the existence of superconducting transition in the two-dimensional 
Hubbard model remains unresolved. Pair correlation functions had been calculated 
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by using QMC methods, and their results were negative for the existence of 
superconductivity in many works. The objective of this paper was to reexamine this 
question by elaborating a sampling method of quantum Monte Carlo method. 

We have calculated the d-wave pair correlation function Pd for the 2D Hubbard 
model by using the QMC method. In the half-filled case Pd is suppressed for the 
repulsive U > 0, and when doped away from half-filling iVe < N, Pd is enhanced 
slightly ioT U > 0. It is noteworthy that the correlation function Pd is indeed 
enhanced and is increased as the system size increases in the 2D Hubbard model. 
The enhancement ratio increases almost linearly oc L as the system size is increased, 
which an indicative of the existence of superconductivity. Our criterion is that 
when the enhancement ratio as a function of the system size L is proportional to 
a certain power of L, superconductivity will be developed. This ratio dependes on 
U and is reduced as U is decreased. The dependence on the band filling shows a 
dome structure as a function of the electron density. In the 10 x 10 system, the 
ratio is greater than 1 in the range 0.3 < n^ < 0.9. This does not immediately 
indicates the existence of superconductivity. The size dependence is important and 
is needed to obtain the doping range where superconductivity exists. Let us also 
mention on superconductivity at half-filling. Our results indicates the absence of 
superconductivity in the half-filling case because there is no enhancement of pair 
correlation functions. 

We have compared two methods: diagonalization QMC and Metropolis QMC. 
For small systems, the results obtained by two methods are quite consistent. When 
the system size is large, Pd(^) is inevitably suppressed and almost vanishes if we use 
the Metropolis QMC method. Pd(^) decreases as the division number m increases in 
this method. We wonder if this excessive suppression of Pdi^) is true. In fact, the 
correlation function Dyy for the ladder Hubbard model obtained by the Metropolis 
QMC also shows a similar behavior when the size is increased, in contrast to enhanced 
Dyy indicated by the density- matrix renormalization (DMRG) method[45]. The 
results by the diagonalization QMC are consistent with those of DMRG [21]. There is 
a possibility that this has some relation with the negative sign. 
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